#!/usr/bin/perl 
#===============================================================================
#
#         FILE: tandem_ms_frag_lib_08_19.pl
#
#        USAGE: ./tandem_ms_frag_lib_08_19.pl  
#
#  DESCRIPTION: 
#
#      OPTIONS: ---
# REQUIREMENTS: ---
#         BUGS: ---
#        NOTES: ---
#       AUTHOR: Han Hu (hh), hh1985@bu.edu
#      COMPANY: Bioinformatics Program & Center for Biomedical Mass Spectrometry, Boston University
#      VERSION: 1.0
#      CREATED: 8/27/2011 12:03:19 PM
#     REVISION: ---
#===============================================================================

use strict;
use warnings;

use Class::Struct;
use List::Util qw[sum];
use Getopt::Std;

## The module I created for storing constant parameters.
use MassData qw(%GAGInfo %BasicInfo);
## The input should be $outfile, $dirname, $ppm, $pre, $z
my %Options;
getopt('dpozxi',\%Options);
our ($dirname, $ppm, $outfile, $format, $pre_mz, $pre_z);
#my @time_stamp = localtime(time);
#$time_stamp[5] += 1900;
#my $time_string = join "_", @time_stamp[4,3,5];
#
if(!defined($Options{x})){
	die "$pre_mz is undef: The value of -x is compulsory!\n";
} else {
	$pre_mz = $Options{x};
}
if(!defined($Options{d})) {
	$dirname = "./output";
} else {
	$dirname = $Options{d};
}
if(!defined($Options{p})) {
	$ppm = 5e-6;
} else {
	$ppm = $Options{p};
}
## Reformat the output of precursor m/z.
my $pre = int($pre_mz + 0.5);
if(!defined($Options{o})) {
	$outfile = join "_", ("gag_library", $pre);
	$outfile = $outfile . ".csv";
} else {
	$outfile = $Options{o};
}
if(!defined($Options{z})) {
	die "$pre_z is undefined: the value of -z is compulsory!\n";
} else {
	$pre_z = $Options{z};
}

##############################
## 0. Define data structure ##
##############################
## 0.1 Basic information

our @links = @{$GAGInfo{'HS'}->{'link'}};                             # A

### Calculate the mass of structural units
## HexA
our @hexa_group = @{$GAGInfo{'HS'}->{'modification'}->[0]};
## GlcN
our @glcn_group = @{$GAGInfo{'HS'}->{'modification'}->[1]};

## 0.2 Monosaccharide unit and fragments
struct MonoUnit => {
	idx => '$', # Position of the monosaccharide unit.
	mass_piece => '@',
	modification => '$',
	fragment => 'UnitFrag',
};
struct UnitFrag => {
	piece => '@',
	cleavage => 'UnitCleavage',
};
struct UnitCleavage => {
	begin => '$',
	end => '$',
	direction => '$',
};
## Chain id => Fragment id => Unit id.

#################################
## 1. Precursor -> Composition ##
#################################

use PrecursorToComposition;

## Calculate the composition.
my $compo_status = PrecursorToComposition::PrecursorToComposition($pre_mz,$pre_z,$ppm);
our @composition;
if($compo_status == 0){
	die "The composition could not be calculated, please recheck PrecursorToComposition.pm!\n";
} else {
	@composition = @$compo_status;
}

our $dp = $composition[0] + $composition[1] + $composition[2];
print "The composition is: @composition\n";

#################################
## 2. Composition -> Sequences ##
#################################
## Update candidate chain with possible positions
use CompositionToSequence;
## Depth first.
## Update the member 'modification'
my $candidate_chain = CompositionToSequence::CompositionToSequence(\@composition);

#===  CLASS METHOD  ============================================================
#        CLASS: BioChemFilter
#       METHOD: filter
#   PARAMETERS: $BioSyncSwitch
#      RETURNS: ????
#  DESCRIPTION: TBD
#       THROWS: no exceptions
#     COMMENTS: none
#     SEE ALSO: n/a
#===============================================================================

use BioChemFilter;
my $BioSyncSwitch = 1;
$candidate_chain = BioChemFilter::BioSync($candidate_chain);

##############################
## 3. Sequence -> fragments ##
##############################
use SequenceToFragments;

## This number can be altered.
## Further job is to enclose this section.
## Initialization of the possible chain.
my $chain_id = 0;

#-------------------------------------------------------------------------------
#  Internal cleavage still needs to be considered, but only focus on the mass loss instead of the exact position.
#-------------------------------------------------------------------------------
#my $fullfilename = $dirname . . $outfile;

open(FG, ">$dirname/$outfile") or die "Cannot write to file $outfile:$!\n";
#open(LOG, ">temp.log");
my %MassLoss = %{$BasicInfo{'MassLoss'}};
my @neutral_loss = ($MassLoss{'h2o'}, $MassLoss{'sg'}, $MassLoss{'sg'}+$MassLoss{'h2o'}, $MassLoss{'co2'}, $MassLoss{'co2'}+$MassLoss{'h2o'}, 2*$MassLoss{'sg'}+$MassLoss{'h2o'},$MassLoss{'co2'} + $MassLoss{'sg'},$MassLoss{'h2o'} + $MassLoss{'co2'} + $MassLoss{'sg'}, 2*$MassLoss{'sg'});

foreach my $chain (@{$candidate_chain}){

	my $chain_pattern = join '-', map {$_->modification} @$chain;

	my $truncated_frags = SequenceToFragments::EdgeCut($chain,$dp);
	
	## Calculate possible mass loss instead of exact position of cleavage.
	foreach my $seg_frag (@$truncated_frags) {
		## Select qualified id of the fragment units.
		## Terminal should not be taken into account.
		my $internal_loss = SequenceToFragments::InternalCut($seg_frag);

#-------------------------------------------------------------------------------
#  Calculte the number of sulfate groups which is used to decide the neutral loss.
#-------------------------------------------------------------------------------
		my $sulfate_count = 0;
		my @modis = map {$_->modification} @$seg_frag;
		foreach my $i (@modis){
			my @site_modi = split '', $i;
			foreach my $j (@site_modi){
				$sulfate_count++ if $j!=0 && $j!=1;
			}
		}

#-------------------------------------------------------------------------------
#  Prepare possible neutral loss.
#-------------------------------------------------------------------------------
		my $neutral_loss;
		if($sulfate_count==0) {		# No sulfate group
			$neutral_loss = join '-', @neutral_loss[0,3,4];
		} elsif($sulfate_count==1) {
			$neutral_loss = join '-', @neutral_loss[0..4,6,7,8];
		} else {
			$neutral_loss = join '-', @neutral_loss;
		}

#-------------------------------------------------------------------------------
#  Format the output.
#-------------------------------------------------------------------------------
		my $start_cleavage = join '-', ($seg_frag->[0]->fragment->cleavage->begin, $seg_frag->[0]->fragment->cleavage->end);
		my $end_cleavage = join '-', ($seg_frag->[$#{$seg_frag}]->fragment->cleavage->begin, $seg_frag->[$#{$seg_frag}]->fragment->cleavage->end);
		my $start_unit = $seg_frag->[0]->idx;
		my $end_unit = $seg_frag->[$#{$seg_frag}]->idx;
#		foreach(@$seg_frag){
#			for(my $i=0; $i<@{$_->fragment->piece};$i++) {
#				if(!defined($_->fragment->piece->[$i])){
#					die "Stop here: $chain_id\t$chain_pattern\t$start_cleavage\t$end_cleavage\t$start_unit\t$end_unit\n";
#				}
#			}
#		}
		my $mass = sum(map {sum(@{$_->fragment->piece})} @{$seg_frag});

		
		print FG "$chain_id\t$chain_pattern\t$mass\t$start_unit\t$end_unit\t$start_cleavage\t$end_cleavage\t$neutral_loss\t$internal_loss\n";

		#print "$chain_id\t$chain_pattern\t$mass\t$start_unit\t$end_unit\t$start_cleavage\t$end_cleavage\t$neutral_loss\t$internal_loss\n";
	}

	$chain_id++;
	next;
}	

close(FG);
#close(LOG);
